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ABSTRACT 

The fraction of ionizing photons which escape their host galaxy and so are able to ionize 
hydrogen in the inter-galactic medium (IGM) is a critical parameter in studies of the 
reionization era and early galaxy formation. In this paper we combine observations 
of Lya absorption towards high redshift quasars with the measured UV luminosity 
function of high redshift galaxies to constrain the escape fraction (/esc) of ionizing 
photons from galaxies at z ~ 5.5 — 6. We employ an A^-body simulation to describe 
the density and peculiar velocity fields and to identify virialized halos. To model the 
ionizing background (IBG) we associate these halos with star-bursting galaxies using 
a semi-analytic prescription. We extract ensembles of mock Lya absorption spectra 
for a range of values of fcsc assuming all resolved (M > 10^ Mq ) simulation halos are 
able to contribute to the IBG. The observed Lya transmission constrains the escape 
fraction to lie in the range /esc ~ 10 — 25% (at z ^ 5.5 — 6). Excluding halos with 
M < IO^^Mq (as might be expected if galaxy formation is suppressed due to the 
reionization of the IGM) implies a larger escape fraction of /esc ~ 20 — 45%. Using 
the numerical results to calibrate an analytic relation between the escape fraction and 
minimum galaxy halo mass we also extrapolate our results to a mass {M ~ 10^ Mq ) 
corresponding to the hydrogen cooling threshold. In this case we find /esc ~ 5 — 10%, 
consistent with observed estimates at lower redshift. We find that the escape fraction 
of high redshift galaxies must be greater than 5% irrespepctive of galaxy mass. Based 
on these results we use a semi-analytic description to model the reionization history of 
the IGM, assuming ionizing sources with escape fractions suggested by our numerical 
simulations. We find that the IBG observed at z ~ 5.5 — 6 implies a sufficient number 
of ionizing photons to have reionized the Universe by z '^ 6. However, if the minimum 
mass for star-formation were > 10^ Mq , the IBG would be over-produced at z < 5. 
In summary, our results support a scenario in which the IGM was reionized by low 
mass galaxies. 

Key words: cosmology: diffuse radiation, large scale structure, theory - galaxies: 
high redshift, inter-galactic medium 



1 INTRODUCTION 

The sources thought to have produced the UV radia- 
tion which reionized the hydrogen gas in the inter-galactic 
medium (IGM) a r e star -bursting galaxies and quasars [e.g. 
iBarkana fc Loebl (|200ir )]. However, the number density of 
the quasar population is observed to decline exponen- 
tially with redshift at z > 2.5, implying that galax- 
ies probably contribute the bulk of UV photons which 
drive reionization at z > 6 and also the evolution of 
the ionizing background post-overlap between z ~ 3 — 6 



Madau et al.' '1999; 'Fan et alj |2002| : ISrbinovskv fc Wvithd 
2007 : Bolton fc Haehnelt 20071 '). The potential contribution 
of galaxies to the UV radiation field is a function of their 
number density, the star formation rate (SFR) within the 
galaxies and the spectral energy distribution (SED) of the 
stellar population. However, the fraction of ionizing photons 
which escape their host galaxy into the IGM, is limited by 
intervening absorption in the inter-stellar medium (ISM). 
The value of this escape fraction (/osc) is therefore a critical 
parameter in studies of reionization, and remains relatively 
uncertain. 

Owing to its importance, constraints on the value of the 
escape fraction feature extensively in the literature. How- 
ever, the definition of fcsc varies depending on the study. 
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A tendency exists in observational studies to quote values 
of /esc as the escape fraction of Lyman- limit photons, rela- 
tive to the escape fraction at 1500A (fcsc")- In this paper, 
we define /esc as the absolute escape fraction of photons 
at the Lyman- limit, since this is the parameter of physical 
importance for the reionization history of hydrogen. Where 
necessary, we convert values of the escape fraction from pub- 
lished studies to be consistent with this definition, assuming 
fcsc" ~ 0.24 and an extinction with a median value of E(B- 
V) = 0.15 at 2 < 3 (|Sianaet al.ll2007l ). 

Early theoretical studies predicted values for the escape 
fraction of between U^ ^ 3 - 15% llDove fc ShuJil 1 1994 



iDove et al.ll2000l:ICiardi et al.ll2002l : lFuiita et al.ll2003l ). how 



■ Clarke fc OevI l|2002l) suggested that /esc could be much 
higher than this if the feedback de pendent poro s ity of the 
ISM is considered. At high redshift, iFuiita et all (|2003l ) ar- 
gue that increased supernovae activity create escape tunnels 
through the ISM for ionizing radiation, and estimate /esc > 
20% by z > 5. Similarly, iRazoumov fc Sommer-LarsenI 
(|2006t ) find that /esc evolves from ~ 1 - 2% at z = 2.39 
to ~ 6 - 10% at 2 = 3.6. In contrast, IWood fc Loebl (|200G| ') 
suggest that /esc decreases with redshift (becoming ~ 1% at 
z ~ 10) due to the inc r easin g den sity of the galactic d isk. At 
z ~ 20.'Wha lenet al.l l|2004l ') and lAlvarez et al.1 (120061 ) . have 
derived escape fractions as high as /esc ~ 95%, however in 
these cases the large value can be attribute d to the assumed 
popul ation-Ill stellar population. Finally, iKitavama et al.l 
(I2OOJ) studied the evolution of Hii regions surrounding a 
central population-Ill star embedded in a galaxy and find 
that /esc increases with the mass of the star and inversely 
with the mass of the halo. 

In the local Universe, observational estimates of 
/esc ~ 2—11%, have been deriveqj from measure- 
ments of e mission line str e ngths r elative to the Lyman - 
limit flux iLeitherer et al.1 ll 19951): iHurwitz et all lll997h: 



iHeckman et al.l (|200ll ): iBergvall et all l|2006h : ICrimes et all 
1 20071) ]. An alternative approach measures the flux in the 
Lyman-continuum relative to the flux at a more readily de- 
tectable frequency above the Lyman-limit (usually 1500A). 
Comparison is then made to the intrinic flux ratio predicted 
by an assumed SED of galaxies, which can be characterized 
by the flux decrement acr oss the Lyman break (LB). As- 
suming LB = 3, Dcharvc ng et al.l ([2( 301) estimates values 
of /esc ^ 10% in nearby galaxies, and iMalkan et al.l (|2003l ) 
estimates values of /esc ^ 1 — 6% from a sample of 11 blue 
galaxies at z ~ 1. At z ~ 3, ISteidel et al.l l|200ll ) compile a 
composite spectrum from 29 Lyman-break galaxies and find 
an average escape fraction of /esc > 10%, and from observa- 
tions of 2 galaxies ICiallongo et al.l l|2002l ) determine values 
of .fesc 5, 4%. From a sample of 14 star-forming galaxies, 
IShaplev et al.l (|2006l ) determines the escape fraction to be 
/esc ~ 3%, however, in two of these galaxies detection of Ly- 
man continuum flux suggests an escape fraction of the order 
of 10% . A less stringent constraint is found bv llnoue et al.l 
l|2005l) . who flnd the escape fraction to be /esc ;$ 20 — 40% at 
« ~ 3 using narrow band images of two galaxi es. Fina l ly, the 
escape fraction has also been constrained bv ISchavd (|2006l ) 



^ In this range fo r /esc we have excluded the finding of 
iHurwitz et al.l Ill997|) for one of the four galaxies in their stud- 
ied sample in which /esc ^ 57%. 



who found /esc ~ 10% at z ^ 3, by comparing the ab- 
sorption systems near an ionizing source with those in the 
average IGM. 

Alternatively, if an SED is available for the galaxy, then 
the the effects of intrinsic and intervening absorption can be 
modelled. Using this method iFernandez-Soto et al.l (|2003r ) 
find /esc ^ 4%, by fitting to the obser ved SEDs of 27 gala xies 
at 1.9 < 2 ^ 3.5. Similarly, at 2 ~ l. lSiana etal] (|2003) fits 
optical/near-infrared model SEDs to the observed spectra 
of 21 galaxies and finds the average, intrinsic Lyman-break 
factor to be ~ 6, from which they conc lude the global es - 
cape fraction to be, /esc ^ 4%. Recentlv. IChen et al.l (|2007l ) 
estimate /esc ^ 2% &t z ^ 2 from the observed distribu- 
tion of neutral hydrogen column densities in the after-glow 
spectra of long duration GRBs. This approach alleviates the 
dependence on largely unknown intrinsic spectral properties 
of the galaxy, however presents another issue of a possible 
GRB environment bias. Finally, combining estimates of /esc 
from published studies with values derived fro m measure- 
ments of the ionizing background, llnoue et al.l (2006) find 
that the escape fraction increases from /esc ~ 1 — 10% in 
the redshift interval 1 < z < 4. 

While there is as yet no consensus on the value of /esc, 
the implications for the reionization history have been ad- 
dressed b y several authors. R ecently, detailed numerical sim- 
ulations i Gnedin et al.ll2007l ) have predicted a value for /esc 
of between 1 and 3%, from halos of mass M > 5 x lO^^M© 
within the redshift range 3 < 2: < 9. In addition to a s mall 
escape fraction in massive galaxies I Gnedin et al.l (|2007|) fur- 
ther predict that halos with Af < 5 x 10^" M© have an escape 
fraction which is negligibly small. This very low efficiency of 
reionization would have profound effects on the reionization 
history. Indeed,, Gnedin (2007) argues that the small escape 
fraction for lower mass galaxies implies that the observed 
galaxy population at 2 ~ 6 could not have produced enough 
UV photons to reionize the IGM. This suggests that the 
process of reionization takes place in a "photon starved" en- 
vironment. This co n clusio n has been previously reached by 
iBolton fc Haehneltl (|2007l ) who studied the ionizing back- 
ground (IBG) based on the Lya absor ption properties to- 
wards high redshift quasars. However, iBolton fc Haehneltl 
(|2007|) inferred a higher a value for the escape fraction of 
/esc ~ 20 -30% (at 2~6). 

In this paper we emp loy an approach that is similar to 
IBolton fc Haehneltl (|2007l ). and constrain the escape frac- 
tion of ionizing photons (at z ~ 5.5 — 6) by combining the 
observed luminosity function (LF) of galaxies with the ob- 
served opacity of the high redshift IGM. However, with re- 
spect to the study of the escape fraction our modelling in- 
cludes several improvements. We employ an A''-body simula- 
tion to describe the density and peculiar velocity fields and 
to identify virialized halos. We determine the UV luminosity 
of each halo using a semi-analytic description and constrain 
model parameters by fitting to the observed galaxy LF. We 
then generate an inhomogeneous ionizing background using 
these discrete galaxies, which depends on the escape frac- 
tion of ionizing radiation. For a range of escape fractions we 
extract an ensemble of mock Lya absorption spectra along 
3072 sight-lines through our simulation box and determine 
the value of /esc which best models the observed Lya trans- 
mission. We note that in our model the escape fraction is not 
explicitly dependent on galaxy mass or angular direction. 
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The paper is organized as follows. In § [2] we outline 
our models for the IBG due to discrete galaxies. We also 
describe our computation of the Lya photon transmission. 
In § [3] we show our results for /eac at several redshifts as- 
suming MFPs which are consistent with the observation- 
ally derived values from Fan et al. (2006). In § 0] we re- 
peat our analysis to constrain /esc assuming the suppres- 
sion of galaxy formation in low mass halos (M < IO^^Mq), 
while in § [S] we exclude halos with M < 5 x 10^*^ Mq. 
In § [6] we extrapolate our numerical results to other min- 
mum galaxy masses. In § [7] we discuss the implications 
of our results using a semi-analytic model of reionization 
history. Finally we make some concluding remarks in § [8] 
Throughout the paper we assume a ACDM cosmology with 
parameters (erg, /i, Ob, ^m, fiA) = (0.8, 0.7, 0.04, 0.3, 0.7). All 
distances are in co-moving units unless otherwise stated. 



2 MODELLING Lya TRANSMISSION AND A 
FLUCTUATING IONIZING BACKGROUND 

A brief summary of our approach to constrain the escape 
fraction is outlined below, followed by several subsections 
describing our model in more detail. We employ an A'^-body 
simulation to compute the density and peculiar velocity 
fields and also to locate dark matter (DM) halos, which we 
then associate with galaxies. To a fraction of these galaxies 
which are considered to be actively star-bursting, we assign a 
luminosity (at A = 1350A) based on models for the SFR and 
the intrinsic SED. We constrain parameters in these models 
by measuring the LF from the simulation and fitting to the 
observed LF. Assuming the IBG to be dominated by galaxies 
(see § H]), we next compute the luminosity of these actively 
UV luminous galaxies at the Lyman-limit, and using a sim- 
ple model generate a field of ionizing radiation which reflects 
the distribution of these discrete sources. We are then able 
to compute an ensemble of Lya forest spectra along lines- 
of-sight through our simulation volume, from which we can 
measure the mean transmission of the IGM. Our model for 
describing the IBG contains the parameter /osc, which we 
adjust so that the mean transmission measured from our 
model matches the observed mean transmission. 



2.1 The A'^-body simulation 

The simulation followed 1024"^ DM particles in a cubic vol- 
ume of side-length, L = 65.6 Mpc//i. Snap-shots of the den- 
sity and peculiar velocity fields were taken periodically in 
redshift space with bound halos identifled using a friends-of- 
friends algorithm. The simulation accurately resolves halos 
with mass, M > 2 x IO^Mq. We place the DM halos and the 
density and velocity fields onto a 512"^ grid. In this paper we 
consider redshift snap-shots of the density and peculiar ve- 
locity fields at redshifts z = 5.5, 5.7, 6.0. The numerical sim- 
ulation assume d crs = 0.9 which we translate to a lower value 
of (78 following jMcQuinn et al.l ([2007). F urther detai l s of th e 
numerical simu l ation can be found in IZahn et al.l (|2007|); 
iMcQuinn et~aLl (|2007l ) and iLidz et al] l|2007h . 

Baryonic gas pressure is not included in our (DM only) 
A^-body simulation. However Hydro-Particle-Mesh (HPM) 
simulatio ns can be used to a pproximate the effect of gas 
pressure. iGnedin fc Huil l|l99a) show that HPM simulations 



of the IGM are capable of reproducing measurable quantities 
(such as absorption spectra) with a precision comparable to 
full hydro-dynami c simulations. Thoug h less reliable than 
HPM simulations, iGnedin fc Huil ([1993) also show that the 
effect of pressure can be incorporated in A'^-body simulations 
by smooth ing the linear den sity fleld on a suitable filtering 
scale (fcfl. lLidz et al.l (|2006l ) measured the power spectrum 
of Lya transmission fluctuations [PF(fc)], using HPM simu- 
lations to describe the density and velocity fields. To deter- 
mine the correct filtering scale for our A^-body simulation we 
therefore measure Pf (fc) following smoothing of the density 
field on a range of fil t ering scales. We find the best agree- 
ment with lLidz et al.l l|2006l) is obtained for a filtering scale 
of fc/ = 20 Mpc~^, and adopt this value throughout the 
paper. 



2.2 The ionizing background 

The fraction of galaxies identified in our simulation which 
are actively star-bursting at a particular redshift, can be 
described by the duty-cycle of the star- burst (tDc = is/in, 
where in is the Hubble time at z). To this fraction of galaxies 
we assign a luminosity (a t A = 1350A) by fi rst estimating 
the SFR which [following iLoeb et al.l ( 20051 )] is dependent 
on the mass of the halo (M), the star- formation efficiency 
(/*) and again the duration of the star- burst (ig). 



SFR =0.17 



0.1 



M 



U 



109 A/0 ) V 10** yr 



0.17 



Mq yr \ 



(1) 



The specific luminosity at wavelength A of each halo can 
then be computed from 



i.(A)=JL(A) 



SFR 



1 M, 



O yr' 



— ergs s 



■^Hz- 



(2) 



where Jl(A) is the specific rest frame luminosity, per Hz, 
per SFR (solar-masses per year). We compute Jl(A) using a 
m odel for SED of star-b ursting galaxies which is presented 
in lLeitherer et al.l l| 19991). We a ssum e metal enriched stars of 
0.05 solar metalhcity, a iScald (|l998l ) IMF (for 1 - lOOM© ) 
and continuous star-formation over a star-burst lifetime of 
100 Myr. This low solar metalhcity is reasonable at high red- 
shift and this star-burst lifetime is consistent with the duty- 
cycle which is constrained below. The relative flux density 
between 1500A and 900A for this SED is given by LB ~ 3. 
As our aim is to reproduce the observed LF we correct 
this intrinsic source luminosity by including extinction by 
dust within the galaxy. We adopt a factor of 1.5 for the ra- 
tio of th e intrinsic to the obs erved luminosity (at z ~ 6) 
following iBouwens et al.l (|2007l ) [and references therein]. For 
10, 000 combinations of to and /* we compute the LF in 
units of the number of galaxies per AB-magnitude per co- 
moving volume (dn/dAf 1350 ) and constrain these parameters 
using a maximum likeli hood technique , fitting our model LF 
to the observed LF of iBouwens et al.l (|2007l ') at 2 = 6. We 
assume flat prior probability distributions for the parame- 
ters t-Dc and /*. In Figure [T] we show the best fit LF. For 
display purposes a line joins the model points (solid) which 
are the centers of the constructed luminosity bins. At the 
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faint end of the LF we have constructed an additional bin 
as within this star-formation model, the simulation resolves 
galaxies down to a lower luminosity than are currently ob- 
served. The error bars on the model reflect the number of 
galaxies in each luminosity bin. The open c ircles represent 
the d ata from the observed UV galaxy LF (|Bouwens et al.l 
[20071). The model points are shifted slightly to the right of 
the observational data for clarity. Inset in Figure [l] is a plot 
of contours representing the 1, 2, 3 — o" confidence levels on 
the parameters /* and ioc- A cross marks the point of maxi- 
mum likelihood, with values of tDc = 0.16 (ts ~ 2.2 x 10* yr) 
and /* = 0.11. 

In this model the minimum mass halo {M ~ 2 x 
IO^Mq ) which contributes to the IBG corresponds to a spe- 
cific luminosity of L ~ 8.9 x 10^*^ ergs/sec/Hz (at 1350 A) 
and a magnitude of M1350 ~ —15.8. The largest halo identi- 
fied in our simulation has a mass of M ~ 3 x 10^^ Mq , which 
given the constrained parameters corresponds to a specific 
luminosity of L ~ 1.2 x 10'^" ergs/sec/Hz and a magnitude 
of M1350 ~ —23.8. However at z = 6, only ~ 13% of iden- 
tified halos have masses greater than 1 x W^'^Mq and only 
~ 1% have masses greater than 5 x 10^° Mq . The largest 
actively star-bursting halo selected in our simulation has a 
mass of M ~ 1 x IO^^Mq , which corresponds to a spe- 
cific luminosity of L ~ 4.4 x W^^ ergs/sec/Hz and a mag- 
nitude of M1350 ~ —22.5. We note that in our model the 
minimum luminosity halo that contributes to the UV back- 
ground is equivalent to the definition of the minimum lumi- 
nosity cut off in the integration of the LF which is used to 
deri ve the global emissivity du e to galaxies in other studies 
[e.g. iBolton fc Haehneltl l|2007l )]. 

Using the constrained values of /* and toe, we re- 
evaluate the luminosity of each star-bursting galaxy at the 
Lyman-limit (912^4). We then re-position the UV luminous 
halos at the center of the simulation cell in which they reside 
and construct a luminosity grid [i(x)] as a function of posi- 
tion, which matches the dimensions of our density and veloc- 
ity grids from the simulations. As the grid is relatively fine 
(65.6/512 Mpc/h) the effect of this re-positioning on the re- 
sulting IBG is negligible in comparison to other uncertainties 
in our model (such as the assumed MFP) . More significantly 
however, the sources remain correlated with the over-density 
field, which is ultimately the desired effect here. Very occa- 
sionally more than one halo resides within the boundaries 
of a cell. In this case we consider a single more massive halo 
at the center of the cell. 

Given the positions and luminosities of the galaxies we 
are able to compute the flux at every positio n on our grid 
using a procedure which is simila r to that in Bolton et al.l 
l|2005h : [Boiton fc Haehneltl (|2007^ : ^McDonald et al.l (|2005h . 



We flrst compute the co- moving energy density (^) at the 
Lyman-limit through the Fourier convolution 



dV 



(x) 



1 



dx L{x.)G{yi' - x) 



(3) 



which conveniently accounts for flux contributions from 
galaxies across the periodic boundaries in our simulation 
box. In this equation Vcdi is the co-moving volume of a grid 
cell and the attenuation of photons as they propagate from 
sources is described by the filter function G(x' — x), which 
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Figure 1. The best fit model LF with parameter values of 
tnc = 0.16 (ta = 2.2 X 10* yr) and /* = 0.11, as constrained 
by a maximum likelihood analysis. Inset is is a plot of the likeli- 
hood distribution in which the contours represent the 1, 2,3 — tr 
confidence levels with a cross to mark the point of maximum 
likelihood. For display purposes a line joins the model points 
(solid). The error bars on the model reflect the number of galax- 
ies in each luminosity bin as measured from our simulation. The 
open circles represent the data from the observed UV galaxy LF 
IjBouwens et al.l 12007! ). The model points are shifted slightly to 
the right of the observational data for clarity only. 



is defined as 



G'(x' - x) 



47r(x' — x)^c 



(4) 



where A„,ip is the MFP of ionizing photons, x' is the position 
at which the energy density is to be evaluated, x is the posi- 
tion of the source, and c is the speed of light. We then com- 
pute the flux in physical units of J21 (ergs/Hz/s/cm^/sr) 
from the co-moving energy density. 



J2i(x ) = -— ^ 7-:H7(x ) 



10- 



A-K dV 



(5) 



Here, the (free) parameter /esc represents the fraction of ion- 
izing photons produced during star-formation which escape 
the host galaxy. An important simplification in the construc- 
tion of our IBG is the assumption of a global, isotropic es- 
cape fraction and a universal MFP at each redshift. 

At z = 5.5, 5.7 , 6.0, w e take values of the MFP esti- 
mated by iFan et al.l (|2006l ) , which assume a frequency av- 
eraged ionization cross-section, (cr^), corresponding to an 
IBG spectrum defined by J^ oc v~^ . We re-compute (a^) us- 
ing the spectrum appropriate for our model and adjust the 
values for A^fp accordingly. At the above redshifts we con- 
sider the MFPs, A„tp = 28.7, 18.9, 4.2 Mpc//i respectively. 
Although these MFPs are derived from observational data, 
we note that they are model dependent and assume a uni- 
form IBG. We discuss the uncertainty in the MFP further 
in? 13:21 
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2.3 Lya absorption and the ionizing background 

Our simulation specifies (as a function of position) the den- 
sity and tiie peculiar velocity of the gas as well as the ion- 
izing radiation field. We c onstr uct mock Lya forest spectra 
as described in lHui et al.1 (|l997l ) and briefiy outline only the 
main points here, referring the reader to the original paper 
for details. Along a gi ven sight-line, over a distance between 
x^and Xr. IHuI et al.l (|l997ii) compute the optical depth as a 
function of observed frequency [t{vo)] as 



limit given by 



^{^o) = 



dx 

TTi 



nHi(x)CTo 



(6) 



Here z is the redshift at the center of the 1-D spatial interval 
Xb — Xa in which dx is the distance between cells, aa is the 
Lya cross-section and nHi(x) is the number density of neu- 
tral hydrogen. We consider sight-lines parallel to the edges 
of the simulation box in which the perpendicular, central 
plane of the box corresponds to 2, which is also the redshift 
of our simulation snap-shot. 

The LyQ cross-section is expressed as a function of the 
velocities uo and u which can be interpreted in terms of the 
observed, the intrinsic and the resonant Lya frequencies 






(7) 



where <Jc.,o = 4.5 x 10"^'* cm^ and c is the speed of light. The 
parameter h = {2kBT /mpY'^ is a function of the gas tem- 
perature (T) which accounts for the thermal broadening of 
the line profile. Here ks is the Boltzm ann constant and nip is 
the mass of a proton. iHui et al.l l| 19971 ) find that the IGM gas 
temperature can be described locally by T = To(l + 5b )~'~^ , 
with values of 1.2 < 7 < 1.7. However. iLai et al.l l|2006l ') find 
the impact of temperature fiuctuations on Lya forest spec- 
tra is small, therefore we assume To — 1.7 x IG^K and an 
isothermal IGM (7 = 1). We note that the uncertainty in- 
troduced by neglecting temperature fluctuations is absorbed 
into the uncertainty of our estimated value for the required 
ionization rate. Along a sight-line, uq is determined in each 
pixel by the redshift relative to 2 and the corresponding ve- 
locity acquired due to the Hubble flow. The parameter u 
is determined by the peculiar velocity and the redshifted 
velocity of the Lya resonance. 

As a function of position we compute the proper number 
density of neutral hydrogen as 



nHi(x) = /„Hi(x)nH(x) 



(8) 



where fn-mi'x.) is the neutral fraction and the density of 
hydrogen is given by nH(x) = nHA(x). The relative over- 
density [A(x) = /9(x)/p] is obtained from the A'^-body sim- 
ulation and nn is the spatially averaged number density of 
hydrogen. We compute the neutral fraction (following Hui 
et al. 1997) as 



/"HI \^) 



X 10" 



Jhi(x) 

0.5 



f T 



VlO'' K 

A(x) 



H-2 



0.0125 
3 



(9) 



In equation ((9} the quantity Jhi(x) is the average flux of ion- 
izing photons in the IBG at wavelengths below the Lyman- 



f°°4nJ,-^di' 
Jh:(x) = ^''- 






-diy 



(10) 



where Ji, is the specific intensity as a function of frequency 
and ai, is the Hi photo-ionization cross section. For the spec- 
trum of star-burst galaxies employed in this paper (§ 12. 2p . 
the ionizing fiux can be expressed relative to the specific 
intensity at the Lyman-limit as Jhi — O.7J21 [which we 
computed in equation ([S}]. We note that the numerator in 
equation (|10|) is just the ionization rate [equation (lllf) ]. high- 
lighting the proportionality between the ionizing flux (Jhi), 
the ionization rate (r_i2) and the flux at the Lyman-limit 
(J21) via the SED. 

We note that as it enters into the calculation of the op- 
tical depth [equations (I6I10|| ]. /esc strictly implies the escape 
fraction of all ionizing photons. Recent calcu lations of the 
frequ ency dependence of the escape fraction IjGnedin et al.l 
[2007]) show that /esc increases by a factor of (less than) 2 
between the Hi (13.6 eV) and Hen (54 eV) ionization thresh- 
olds. This suggests that interpreting /esc as the escape frac- 
tion at the Lyman-limit, leads to a small over-estimate of 
the quantity which is constrained via observations of galax- 
ies as summarized in § [l] However, the escape fraction at 
frequencies above the Lyman-limit are weighted by a factor 
proportional to the adopted SED. Therefore, we expect the 
difference between /esc at the Lyman-limit and the escape 
fraction which is explicitly constrained by our model to be 
small in comparison to other uncertainties in our model. 

Finally, the transmission at a particular frequency can 
be related to the optical depth by F{vo) = exp[— r (i/o)]. 
Using the model for computing the optical depth presented 
above, we can construct Lya-forest spectra by computing 
the transmission over a range of {vo) frequencies . 



2.4 The observed transmissivity of the IGM 

For a sample of 19 Lya absorption spect ra observed towards 
high redshift quasars, iFan et al.l (|2006l ) measure the aver- 
age (or effective) transmission over a wavelength range in 
the Lya forest which corresponds to a redshift interval of 
A2 = 0.15. From thes e transmis s ion m easurements along 
the available sight-lines. iFan et al.l (|2006r ) compute the mean 
effective transmission {F°^f') in redshift slices (of thickness 
0.2) of the IGM, which are centered on slightly different 
redshifts to the redshifts of our simulations (2). We there- 
fore estimate values of F"^;" at 2 = 5.5, 5.7 , 6.0 di rectly from 
the measurements of F^a from iFan et al.l l|2006r ) (their Ta- 
ble 2) within redshift bins of 0.2 centered on each 2. We 
compute the mean and the variance in this sample assum- 
ing a Gaussian distribution for uncertainties on individual 
transmission measurements. At 2 = 5.5, 5.7, 6.0 the observed 
mean effective transmission values (and 1 — a scatter) are 
F,°t= = 0.079 ± 0.035, 0.047 ± 0.031, 0.006 ± 0.005. 

We generate an ensemble of 3072 absorption spectra 
for evenly spaced sight-lines throughout our simulation box. 
Along each sight-line we then compute the effective trans- 
mission (-Fcft) over stretches of the forest which enables di- 
rect comparison with the observational data. From the en- 
semble of foH we calculate the mean effective transmission, 
F^is. Finally, we rescale the optical depth in each pixel (via 
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Figure 2. Our estimates of the Hi ionization rates at 2 = 
5.5,5.7,6.0 assuming MFPs A„fp = 28.7,18.9,4.2 Mpc/h and 
mean transmission values of _Fcft = 0.079, 0.047, 0.006, respec- 
tivel y (open circles). For comp arison we show estimates o f F— 12 
from lBolton fc Haehnelt bOOTl ) (diamonds) and lFan etahl (I2OO6I ') 
(points) which are also based on IGM Lya transmission. We also 
show our constrained values of F_i2 assuming uniform IBGs at 
each redshift (open squares). 



/esc) so that our simulated value of F^a matches the observed 
value. 



3 THE IONIZATION RATE AND THE 
ESCAPE FRACTION 

In this section we present our estimates of the Hi ioniza- 
tion rate (F) and the escape fraction of ionizing photons 
(./esc), which we have constrained by adjusting the level of 
the IBG in our model so as to reproduce the transmission 
measurement s (F"^') inferred from the spectra of high red- 
shift quasars (JFan et al.ll2006l) . 



3.1 Constraints on the ionization rate 

The Hi ionization rate is computed from the flux in the IBG 
(Jix) at all frequencies above the Lyman limit (i^l) by 



r — 4-K I -—^a^du 



(11) 



where a^, is the Hi photo-ionization cross-section, and hp 
is Planck's constant. For our spectrum of an IBG due to 
star-burst galaxies this relationship can be expressed as 
r_i2 = 2.8J21, where F_i2 is the photo-ionization rate in 
units of lO^^^s^^. To find the ionization rate in the IGM 
from the Lya forest spectra we trial values of /esc, which is 
the free parameter we use to scale the intensity of the IBG. 
For each /esc, we use equation ((SJ to find J21, and hence 
the neutral fraction via equation ([Qj. We compute Lya for- 
est spectra from the absorption at a range of frequencies 
computed via equation ©. A value for the effective opti- 
cal depth is computed from the Lya forest spectra along 
many sight-lines for a given /esc- The value of /esc is then 



adjusted until the simulation reproduces the observed, mean 
Lya absorption. 

In Figure [2] we show the mean ionization rates con- 
strained using our model for both a non-uniform back- 
ground (open-circles) and a uniform background. We note 
that in our model the ionization rate is spatially non- 
uniform as it follows the non-uniform IBG and so our es- 
timates of the mean ionization rate are explicitly {F_i2(x)). 
We compare our resul ts to t he ionization rates predicted 
by iBolton fc HaehneJtl l|2007l 'l (diamonds) and iFan et all 
(2006) (points). The 1 — a error in the ionization rate re- 
fiects the 1 — a uncertainty in the observed transmission 
at e ach redshift. Our results are in excellent agreement 
with iBolton fc H aehnelt [2003) and a factor of ~ 2 greater 



than those of [Fan et al. (20 0^). W e ref er the reader to the 
discussions in IBolton et al.l (|2005r ) and IBolton fc Haehnelti 
(120071 ) regarding the likely under-estimation of the ioniza- 
tion rates predicted by the fiuc tuating G u nn-Pe terson ap- 
proximation model assumed by iFan et al.l (|2006f ). It is also 
worth noting that our values for r_i2 assuming a uniform 
background are less than those for a non-uniform back- 
ground, in agreement with t he findings of previous authors 
iGnedin fc HamiltonI (l2002f);lMeiksin fc Whitd (|2004l ): ICroftl 
i200i] and lBohon fc Haehnelti (120071 )]. In our model for a 
fluctuating IBG, high ionization rates are biased to over- 
dense regions and the transmissivity of rare voids remains 
largely unaffected by increasing the level of the IBG. In con- 
trast, raising the level of a uniform IBG change s the t rans- 
missivity of all regions [see IBolton fc Haehnelti (|2007l ) and 
references therein]. On the other hand, the contribution of 
low density regions to the aggregate Lya absorption ob- 
served along a LOS is relatively minor, irrespective of the 
IBG. 



3.2 The escape fraction of ionizing photons 

In our model, we adjust the intensity of the IBG (via the 
parameter /esc) in order to reproduce a given effective trans- 
mission (F). Following this constraint of the escape fraction 
which corresponds to a particular F, we calculate the prod- 
uct Cf = /*/csc- Our previous constraint of the SFR via the 
parameters /* and toe (see § 12. 2p describes the distribution 
-TT^^-iTT- From this distribution we then form the distribu- 
tion Stt by marginalizing over ioc- Given '^^ 



df* 



the probability distribution of fc, 
,, ,, , using 



sion (--^^ 



we compute 
for a particular transmis- 



df* 



dP 



dfc, 



Cf 



dP_ 

'df 



dp 



dfc, 



dP Cf 

df* /ese' 



(12) 



The upper panel in Figure |3] shows the probability dis- 
tributions of ^jf^L at 2 = 5.5,5.7,6.0 assuming IBGs 
in which the MFPs of ionizing photons are, A„fp = 
28.7, 18.9, 4.2 Mpc/h. In this figure Cf corresponds to the 
observed mean effective transmission (i.e F — F^it)- Our 
analysis suggests /esc ~ 10 — 25% over the redshift range 
z ~ 5.5 — 6. Theses distributions represent the statistical 
uncertainty in /esc owing to uncertainty in the modelling of 
the LF. 

The lower panel in Figure |3] shows the probability dis- 
tributions of /esc including the additional uncertainty owing 
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Figure 3. The likelihood distribution -rf — . We have adjusted 

"/esc 

the intensity of the IBG such that the model values for the mean 
transmission match the observed values. At z = 5.5,5.7,6.0, 
these mean transmission values are Fdi = 0.079, 0.047, 0.006, 
respectively. The IBG model is defined by assumed MFPs of 
Amtp = 28.7, 18.9, 4.2 Mpc/?t respectively. The upper panel shows 
the distributions of fcac corresponding to the mean transmission. 
The lower panel shows the distributions of /esc after integration 
over the uncertainty in the mean transmission. 



to measurement error in the transmission. To find this dis- 
tribution, 
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Figure 4. Likelihood distributions for /esc and (A/esc) corre- 
sponding to the observed mean transmission at 2 = 5.5 for IBGs 
defined by MFPs 3.5,7.0,18.9,28.7,70.0 Mpc/h. We have ad- 
justed the intensity of the IBG such that the model values for 
the mean transmission match the observed value. At z = 5.5, the 
observed mean transmission is Fdi = 0.079. In the upper panel 
the maximum likelihood occurs for /esc = 70,30,13,11,9% re- 
spectively. In the lower panel the maximum likelihood occurs for 
An,tp/csc = 2.4,2.1,2.5,3.2,6.3 Mpc/h respectively 



dP 

d/esc 



dCl 



dP 

dfcec 



dP dP 

dPdC^ 



(13) 



we assume the observed transmission to follow a Gaussian 
distribution ( ^ ) with values for the mean and variance that 



are given in § 12.41 The distributions. 



dP 

dfcsc 



can be eval- 



uated from equation (|12|) . and -^^ is computed from the 
simulations. The uncertainty in transmission dominates the 
uncertainty in estimates of /esc- We find a statistical relative 
error on the value of /esc that is of order ~ 50%, which is 
larger than the difference between the most likely values at 
each redshift. 

The value of the MFP warrants further discussion. The 
MFPs assumed in our model are based on those derived in 
iFan et al.l (120061 ) which are generally lower (by a facto r of > 
2) than those sugg ested bvlBohon fc Haehneltl (120071 ). The 
MFP estimates of iBolton fc Haehneltl (|2007| ) assume that 
photons propagate freely in regions below some critical over- 
density which may be determined from the ionization rate 
and the col umn density (A'^ij^i ) of self-shielded neutral clouds. 
In contrast. I Fan et al.l (|2006t ) estimate the MFP considering 
the distribution of the neutral fraction over all over-densities 
in the IGM. Althoug h our estimates of the ion ization rate 
are in agreement with iBolton fc Haehneltl (|2007|) , we choose 
to adopt lower values for the MFPs for two reasons. Firstly, 
the MFP of ionizing photons may be over-estimated by a fac- 
tor of ~ 2 by neglecting Hi absorber s between Lyman -limit 
systems (Mhalda- Escude 2003 ; Furlanetto fc Ohl(2005l '). Sec- 
ondly, as the MFP is very uncertain, adopting the lower val- 



ues is conservative with respect to the photon budget for 
reionization as we will discuss in § [T] 

We note that uncertainty in the MFP represents a large 
source of uncertainty in our estimates of the escape frac- 
tion, and we therefore next investigate the sensitivity of the 
results to the MFP. In Figure |4] we plot the distribution 
of the escape fraction (upper panel) corresponding to the 
observed mean transmission at z — 5.5, assuming several 
MFPs (3.5,7.0,18.9,28.7,70.0 Mpc/h). An increased MFP 
leads to a decrease in the escape fraction required to repro- 
duce the same IBG. For the shortest (longest) MFP con- 
sidered [A — 3.5(70.0) Mpc/h] the most likely value for the 
escape fraction is /esc ~ 10.2%(1.3%). This range represents 
a factor of ~ 10 variation in the escape fraction given a range 
of a factor of ~ 20 in MFP. In the lower panel we plot the 



distribution of 



dP 



d(A/cBo) ' 



where 



dP 



dP 



rf(Amfp/esc) df* Amfp/Isc 



(14) 



where Cf corresponds to the mean observed transmission. 
This demonstrates the level of degeneracy of the prod- 
uct A,„fp/csc. For the shortest (longest) MFP considered 
[A„tp = 3.5 (70.0) Mpc/h] the most likely value in the prod- 
uct Amfp/esc is ~ 2.4 (6.3) Mpc/h. This represents a factor 
of less than 3. Thus we can summarize our results with the 



r^ 



where 2 < Kx < 6. 



relation /esc(A„,fp) ~ Kx (m^^/tj^ 

Excluding our estimate of /esc at the longest MFP this range 

for Kx narrows to 2 < Kx < 3. 
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3.3 Further discussion 



Before proceeding it is illustrative to c ompare estimates of 
the e s cape fraction whicli a re derived by (jBolton fc Haehneltl 
bOOTl l- lBolton fc Haehneltl ([2007,) find the escape fraction to 
be 20 — 30% by comparing the constrained ionization rate 
to the global emissivity, which they compute by integrat- 
ing the observed UV luminosity function down to a mini- 
mum galaxy brightness of Afiaso — —18. They also consider 
quasars to contribute to the global emissivity, however this 
contribution is less than ~ 10% and not considered further 
here. In our model we include sources down to a lowe r lu- 
minos ity (Mi35o ~ —16) than used in Bolton & Hae hneltl 
(|2007l ). which corresponds to a factor of ~ 2 difference in the 
integrated emissivity. However this factor i s partially can- 
celled by including the revised galaxy LF of iBouwens et al.l 
l|2007l) , which leads to a decrease in the integrated emis- 
sivity of a factor of ~ 2. In a ddition, the SEP a s sumed 
in our study is ~ 35% harder [Bolton fc Haehneltl ([20071 ') 
use r_i2 — 2.I J21]. Finally, the definiti on of the escape 
fraction used in iBolton fc Haehneltl ([20071 ) is equivalent to 
the relative escape fraction which is frequently quoted in 
observational studies (see § [T]). This is because they infer 
the escape fraction of ionizing photons having assumed an 
emissivity which is normalized by integrating over the ob- 
served LF of LBG galaxies. Following Bouwens et al. (2007) , 
this emissivity should be corrected by a factor of 1.5 to de- 
scribe the intrinsic emission of galaxies (at 1350A) at z ~ 6. 
An important point of d ifference between our work and 
iBolton fc Haehneltl (|2007| ) is that the identification of ha- 
los as the sources allows the escape fraction to be connected 
with halo mass as well as luminosity, rather than luminosity 
alone. 

The above corrections should be considered in order 
to quantitatively c o mpar e our estimate of /esc to that of 
IBolton fc Haehnelt (12007 ). Follow ing these corrections the 
result of lBolton fc Haehneltl l[2007^ becomes /esc ~ 15-25%, 
which is consistent with our estimates. Although both stud- 
ies are normalized to the observed LF, a significant distinc- 
tion in our model is that the IBG is produced by discrete 
sources and so spatially non-uniform. However, th e consis- 
tency of our results with IBolton fc Haehneltl ([2007|) suggests 
that distinction between a global escape fraction derived by 
comparing to the integrated emissivity, and a globally aver- 
aged /esc which is apriori attributed to individual galaxies, 
is not significant. 



Finally, we note an additional source of uncertainty in 
our model which is associated with the assumed isotropic 
esca pe of ionizing rad i ation from galaxies. Previous stud- 
ies jFuiita et al.l l2003l : iRazoumov fc Sommer-LarserJ 120061 : 
iGnedin et al.ll2007l ) have shown that ionizing radiation es- 
capes from host galaxies through ionized tunnels in the ISM 
created by radiative feedback. This may account for the 
range of values derived for /esc from observation of galax- 
ies along different sight-lines. In our model, this anisotropy 
would increase the fluctuations in the IBG. Based on our 
finding that the impact of a non-uniform IBG on the global 
value of /esc is not significant, we expect that the uncer- 
tainty introduced to the value of /esc by the assumption of 
an isotropic escape fraction in our model is not significant. 
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Figure 5. The best fit model LF with parameter values of 
tnc = 0.18 (is = 2.5 X 10* yr) and /* = 0.12, as constrained 
by a maximum likelihood analysis. Inset is is a plot of the likeli- 
hood distribution in which the contours represent the 1, 2,3 — tr 
confidence levels with a cross to mark the point of maximum 
likelihood. For display purposes a line joins the model points 
(solid). The error bars on the model reflect the number of galax- 
ies in each luminosity bin as measured from our simulation. The 
open circles represent the data from the observed UV galaxy LF 
IjBouwens et al.l 120071 ). The model points are shifted slightly to 
the right of the observational data for clarity. 



4 DWARF GALAXY SUPPRESSION 

In the previous section we considered an IBG in which we 
included all halos down to the mass resolution of the sim- 
ulation (2 X 10^ M0). In our best fit model this mass cor- 
responds to a luminosity of M1350 ~ —16. In this section 
we investigate the impact on the inferred escape fraction of 
allowing the suppression of dwarf galaxies. The reionization 
of the IGM elevates the Jeans mass and so increases the re- 
quired minimum virial temperature a DM halo must have 
in order for gas to accrete and form stars. This minimum 
virial temperature could be as large as ~ 2.5 x 10^ K which 
corresponds to a minimum halo mass, M ^ jr ^ 10^° Mq 
jEfstathioul Il992l : iThoul fc Weinberd 1 19961 : iDiikstra et al] 
[2OOJ). The suppression of small galaxies will modify both 
the topology and the level of the IBG and hence also the de- 
rived estimate of the escape fraction. We repeat our analysis 
of the galaxy LF, limiting the mass of halos which are consid- 
ered as potential sources of UV radiation to M > 10^° Mq . 
Repeating our analysis of section 12.21 Figure [5] shows 
the best fit LF for the suppression of galaxies with halo 
mass less than W^^Mq . For display purposes a line joins 
the model points (solid) which are the centers of the con- 
structed luminosity bins. At the faint end of the LF we 
have constructed an additional bin as we are able to re- 
solve galaxies down to a lower luminosity in our simulation 
than are currently observed. The error bars on the model 
reflect the number of galaxies in each luminosity bin. The 
open circles represent the data from the observed UV galaxy 
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Figure 6. The likelihood distribution ft ■ We have adjusted 

"/esc 

the intensity of the IBG such that the model values for the mean 
transmission match the observed values. At z = 5.5, 5.7, 6.0, 
these mean transmission values are F^a = 0.079, 0.047, 0.006, 
respectively. The IBG model is defined by assumed MFPs of 
Amfp = 28.7, 18.9, 4.2 Mpc//i respectively. The upper panel shows 
the distributions of /esc corresponding to the mean transmission. 
The lower panel shows the distributions of /esc after integration 
over the uncertainty in the mean transmission. 
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Figure 7. The likelihood distribution Ji . We have adjusted 
the intensity of the IBG such that the model values for the mean 
transmission match the observed values. At 2 = 5.5, 5.7, 6.0, 
these mean transmission values are F^a = 0.079, 0.047, 0.006, 
respectively. The IBG model is defined by assumed MFPs of 
Amfp = 28.7, 18.9,4.2 Mpc//i respectively. The upper panel shows 
the distributions of /esc corresponding to the mean transmission. 
The lower panel shows the distributions of /esc after integration 
over the uncertainty in the mean transmission. 



LF (JBouwens et al.l [2003). The model points are shifted 
slightly to the right of the observational data for clarity 
only. Inset in Figure [5] is a plot of contours representing 
the 1,2,3 — (T confidence levels on the parameters /* and 
foe. A cross marks the point of maximum likelihood, with 
values of toe = 0.18 {t„ = 2.5 x 10* yr) and /* = 0.12. 
These values are not significantly difi'erent from those pre- 
viously derived for the inclusion of galaxies with halo mass 
down to IO^Mq (§ 12. 2|) . For the suppression of galaxies 
with halo mass less than lO^'^M© , the minimum mass halo 
which contributes to the IBG corresponds to a specific lu- 
minosity of 1/ ~ 4.3 X 10^^ ergs/sec/Hz and a magnitude of 
Mi35o ~ —17.5. This lumin osity is smaller than the mini- 
mum luminosity included in iBolton fc Haehnelti (|2007l ) . 

Finally, we compute the likelihood distributions for the 
escape fraction, analogous to those shown Figure (3) In the 
upper panel of Figure [6] we show the probability distribu- 
tions of /esc at our fiducial redshifts and MFPs, where we 
have constrained the transmission to the mean observed 
transmission. This yields values for /esc between 20% and 
40%, which are larger than those found in Figure [3] where 
we included sources down to a halo mass of ~ 10^ Mq . The 
distributions in the upper panel of Figure [6] represent un- 
certainty in /esc owing to uncertainty in the modelling of 
the LF. The lower panel in Figure [6] shows the probability 
distributions of /esc for the same redshifts and MFPs as the 
lower panel, however here we have included the additional 
uncertainty introduced by the distribution of optical depth 
measurements. As before, we assume the observed transmis- 
sion to follow a Gaussian distribution (^) with values for 
the mean and variance that are given in § 12.41 



5 REIONIZATION BY MASSIVE GALAXIES 



In a recent paper, iGnedin et al.l (|2007|) suggested that the 
escape fraction of ionizing photons is a few percent for mas- 
sive galaxies, but is negligibly small for galajcies which re- 
side in halos of mass M < 5 x 10^° Mq. We investigate the 
implications for the escape fraction inferred from Lya trans- 
mission data by excluding halos of mass M < 5 x 10^^ Mq. 
The suggestion that ionizing photons do not escape galax- 
ies with halo mass less than 5 x 10^" M0 does not imply 
that photons above the Lyman-limit (i.e. 1350A) are sim- 
ilarly affected. Therefore we assume the values for /* and 
trtc from §2] to compute the IBG in this case. In this case, 
the minimum halo mass corresponds to a minimum specific 
luminosity of L ~ 2.2 x lO'^* ergs/sec/Hz and a magnitude of 
A/1350 ~ —19.3, whi ch is brighter than the min imum mag- 
nitude considered bv lBolton fc Haehnelti ([20071). 

Again we compute the likelihood distributions for the 
escape fraction. In the upper panel of Figure [7] we show the 
probability distributions of /esc at our fiducial redshifts and 
MFPs, where we have constrained the transmission to the 
mean observed transmission. This yields values for /esc be- 
tween 46% and 86%, which are larger again than those found 
in Figure [31 where we included sources down to a halo mass 
of ~ 10^ M0 . As previously computed in § 13.21 the lower 
panel in Figure [7] shows the distributions of /esc in which we 
have included the uncertainty introduced by the distribu- 
tion of optical depth measurements at each redshift. These 
results imply that reionization by massive g alaxies alone, 
with escape fractions of order a few percent IjGnedin et al.l 
i2007f ) is not consistent with the measured ionization rate at 
2 ~ 5 - 6. 
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Figure 8. Extrapolation of estimates of fcsc as a function of 
minimum mass from numerical simulations. Dashed lines show 
analytic fits to the numerical data for both the maximum and 
minimum values found for the normalization parameter A (see 
text). The numerical estimates are shown including the range for 
the median value of /osc predicted in our simulation at the three 
redshifts modelled. The analytic fit suggests that the value for fcsc 
is greater than 100% for minimum masses, M^ain ^ lO^'^, and ap- 
proaches a constant value of ~ 5 — 10% for minimum masses near 
the hydrogen cooling threshold. Figure|8]indicates that the escape 
fraction from high redshift galaxies must be > 5% irrespective of 
the minimum mass. 



live variation in this normalization parameter is found to be 
7% (11%) for minimum mass halos of Af = 2x 10^-5x 10^°. 
In Figure[8]we plot /esc(Mmin) (dashed lines) assuming both 
the maximum and minimum values found for A. 

Figure [8] shows that as the minimum mass is increased, 
the global escape fraction increases to compensate for the 
lost emissivity. The analytic fit suggests that the value for 
/esc is greater than 100% for minimum masses, Mmin > 10^^, 
indicating an upper limit to Mmin. The values of /esc(Mmin) 
can be extrapolated to masses below the mass resolution 
limit of the simulation (2 x 10^ Mq ) as shown in Figure [S] 
The escape fraction approaches a constant value of ~ 5—10% 
for minimum masses near the hydrogen cooling threshold. 
These escape fractions are comparable to observed estimates 
for individual galaxies at lower redshift. Interestingly, Fig- 
ure |8] indicates that the escape fraction from high redshift 
galaxies must be > 5% irrespective of the minimum mass. 



7 MODELS OF THE REIONIZATION 
HISTORY 

The Lya absorption probes the instantaneous emissivity of 
the galaxy population, and hence allows the constraints on 
the escape fraction presented thus far. However, the redshift 
at which reionization occurs is dependent on the cumulative 
star-formation. In the remainder of this paper we discuss 
constraints on the reionization history. 



6 EXTRAPOLATION OF NUMERICAL 

RESULTS TO OTHER GALAXY MASSES 

In Figure |8] we summarize our estimates of the escape frac- 
tion for the three minimum mass scenarios we have con- 
sidered. The numerical estimates are shown including the 
range for the median value of fcsc predicted in our simu- 
lation at the three redshifts modelled. For comparison, the 
dashed lines represent simple analytic fits to describe the 
relationship between the minimum mass of halos which con- 
tribute to the IBG and the escape fraction. To estimate this 
relationship we assume that the product of the escape frac- 
tion and the star-formation rate is constant and that the 
star-formation rate can be approximated as being propor- 
tional to the time derivative of the collapsed mass faction, 
dFcoi/dt. The quantity dFcd/dt is computed analytically as 
a f unction of halo mass a nd redshift using the prescription 
of (jSheth fc TormerJll999l ). Thus the escape fraction can be 
expressed as 



fcsc(Mmin) = A{Mrr 



dFcol 



dt 



(Mn 



(15) 



where Mmin is the minimum halo mass of galaxies which 
contribute to the IBG, and A is a constant. 

We normalize the relation between the escape fraction 
and the star-formation rate using our estimates of fcsc from 
the simulations (at z = 5.5, 5.7, 6.0), in which the minimum 
contributing halo masses (M = 2 x 10^, 10^°, 5 x lO"' Mq ) 
are adopted to compute dF/dt. The analytic time deriva- 
tive of the collapsed mass faction [dFcd/dt) does not per- 
fectly describe the SFR from the simulation. The value of 
A is therefore mass and redshift dependent and the rela- 



7.1 Semi-analytic model for reionization 

In this section we use a semi-analytic model to study the 
implications of the implied values of escape fraction on 
the re ionization history of the IGM. iMiralda-Escude et al.l 
(l 200d ) presented a formalism which allows the calculation 
of an effective recombination rate in an inhomogeneous uni- 
verse by assuming a maximum over-density (Ac) penetrated 
by ionizing photons within Hii regions. Their model assumes 
that reionization progresses rapidly through islands of lower 
density prior to the overlap of individual cosmological ion- 
ized regions. Following the overlap epoch, the remaining re- 
gions of high density are gradually ionized. It is therefore 
hypothesized that at any time, regions with gas below some 
critical over-density Ai = Pi/{p) are highly ionized while 
regions of higher density are not. In what follows, we draw 
primarily from their prescription and refer the reader to the 
original paper fo r a detailed discussion of its motivations 
and assumptions. IWyithe fc Loebl l|2003) employed this pre- 
scription within a semi-analytic model of reionization. This 
model was extended by Srbinovsky & Wyithe (2007) and 
by Bolton & Haehnclt (2007). We refer the reader to those 
papers for a full description. 

Within the formalism of Miralda-Escude et al. (2000) 
we describe the post-overlap evolution of the IGM by com- 
puting the evolution of the fraction of mass in regions with 
over-density below Ai, 



FmiAi) 



dAPv(A)A 



(16) 



where Pv(A) is the volume weighted probability distribution 
for A. Miralda-Escude et al. (2000) quote a fitting function 
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Figure 9. Models for the reionization of the IGM and the subsequent post-overlap evolution of the ionizing radiation field. In each 
panel the case shown corresponds to a value for the critical over-density prior to the overlap epoch of Ac = 5. Upper Left Panel: The 
ionization rate as a function of redshift. The observational points are from lBolton &c Haehneltl 1 120071) . Lower Left and Right Panels: The 
volu me and mass averaged frac tions of neutral gas in the universe. The observational points for the volume averaged neutral fraction are 
from IIBolton fc Haehneltll2007l) . while the observed mass-fractions are from the damped Lya measurements of Prochaska et al. (2005). 
Upper Right Panel: The mean-free-path for ionizing photons. The data points are based on Storrie-Lombardi et al. (1994). The models 
are labelled by the assumed values of A/csc and f*^^ . 



which provides a good fit to the volume weighted probability 
distribution for the baryon density in cosmological hydrody- 
namical simulations. This probability distribution remains 
a reasonable description at high redshift when confronted 
with a more modern cosmology and updated simulations, 
although an analytical approximation for the high density 
tail of the distribution remains necess ary as a best guess 
at co rrecting for numerical resolution IjBolton fc Haehneltl 
120071 ). 

Prior to the overlap epoch, a constant value Ai = Ac 
is assumed, while in the post overlap era our model for the 
reionization history computes the evolution of Ai . We com- 
pute evolution of the quantity Qi which is defined to be 
the volume filling factor within which all matter at densi- 
ties below Ac has been ionized. Within our formalism, the 
epoch of overlap is precisely defined as the time when Qi 
reaches unity. However, we have only a single equation to 
describe the evolution of two independent quantities Qi and 
Fm (or equivalently Ac). The relative growth of these de- 
pends on the luminosity function and spatial distribution of 
the sources. We assume Ac to be constant with redshift be- 
fore the overlap epoch and compute results for models with 
values of Ac = 5. Within the formalism of Miralda-Escude et 
al. (2000), a maximum over-density of Ac = 5 corresponds 
approximately to the mean separation between galaxies at 
high redshift. For a large range of values, the effect of vary- 



ing Ac is s maller than other uncertainties in the problem 
JWvithe et al. 2007). 

Our approach is to compute a reionization history given 
a particular value of Ac, combined with assumed values for 
the efficiency of star-formation and the fraction of ionizing 
photons that escape from galaxies. With this history in place 
we then compute the evolution of the background radiation 
field due to these same sources. After the overlap epoch, 
ionizing photons will experience attenuation due to resid- 
ual over-dense pockets of HI gas. We use the description of 
Miralda-Escude et al. (2000) to estimate the ionizing photon 
mean-free-path, and subsequently derive the attenuation of 
ionizing photons. We then compute the fiux at the Lyman 
limit in the IGM due to sources immediate to each epoch, 
in addition to redshifted contributions from earlier epochs. 

As in earlier sections we assume an IBG dominated 
by galaxies and the spectral energy distribution (SED) of 
population-II star forming galaxies, using the model pre- 
sented in Leitherer et al. (1999). The star- formation rate 
per unit volume is computed ba sed on the collapsed frac - 
tion ob tained from the e xtended IPress fc Schechteii (|l974f ) 
model (JBond et al.lll99ll ) in halos above the minimum halo 
mass for star formation, together with an assumed star for- 
mation efficiency. 

In a cold neutral IGM beyond the redshift of reioniza- 
tion, the collapsed fraction should be computed for halos of 



© 2006 RAS, MNRAS 000, [T]-?? 



12 Srbinovsky & Wyithe 



sufficient mass to initiate star formation. The critical virial 
temperature is set by the temperature (Tn ~ 10'* K) above 
which efficient atomic hydrogen cooling promotes star for- 
mation. Following the reionization of a region, the Jeans 
mass in the heated IGM limits accretion to halos above 
7i ~ lO'^ K (Efstathiou 1992; Thoul & Weinberg 1996; Di- 
jkstra et al. 2004b). As described in the introduction, only 
a fraction of ionizing photons produced by stars enter the 
IGM. Therefore an additional factor of /esc (the escape frac- 
tion) must be included when computing the emissivity of 
galaxies. In our fiducial model we assume this escape frac- 
tion to be independent of mass. We define a parameter 



Jci 



ru 



Figure |9] shows our fiducial model for the reionization 
of the IGM and the subsequent post-overlap evolution of 
the ionizing radiation field, assuming /^c = 0.005 (fasc ~ 
5%). In the top left panel of Figure |9] we show the evolution 
of the ionization rate. The observational p oints are from 
the simulations of iBolton fc Haehnelti (120071). which at the 
redshifts shown are based on observations (Songaila 2004; 
iFan et al. 2006) and simulated data (Schaye et al. 2003 ) of 
the Lya opacity of the IGM. The assumed value of /*ac — 
0.005 provides the best fit to the data, and is consistent with 
modelling of the LF (/* ~ 0.1) as well as extrapolation of the 
escape fraction to minimum masses below lO" Mq (/osc ^ 
5%). 

In the lower-left and lower-right panels we plot the 
corresponding volume and mass (upper curves) averaged 
fractions of neutral gas in the universe. The observa- 
tional po ints for the volum e aver aged neutral fraction 
are from iBolton fc Haehnelti l|2007l ). while the observed 
mass-fractions a r e from the damped Lya measurements of 
IProchaska et al.l l|2005l ). and therefore represent lower lim- 
its on the total HI content of the IGM. Both curves for 
the mass fraction and the volume fraction show excellent 
agreement with these observed quantities. This is despite 
their differing by 3 orders of magnitude and indicates the 
applicability of the model over a wide density range. In the 
upper-right panel we plot the evolution of the ionizing pho- 
ton mean-free-path. Th e data points at z > 5 are based on 
IStorrie-Lombardi et al.l (|l994l ). At z > 5 the values of the 
MFP adopted in our simulations are shown for comparison. 
Again the model is in good agreement with the available 
observations and the estima ted MFPs. The observed val - 
ues for the MFP obtained bv lStorrie-Lombardi et all (|l994l ) 
are found from the number density of Lyman-limit systems 
(which represents an upper limit, see discussion in §|3]2]) and 
is independent of the Lya forest absorption derived quanti- 
ties of ionization rate and volume averaged neutral fraction, 
as well as being independent of the HI mass-density mea- 
surements. Our simple model therefore simultaneously re- 
produces the evolution of three independent measured quan- 
tities. 

The agreement of this model with observations is an 
important consistency check of our parameter estimates in 
previous sections. Taken together, our results for the escape 
fraction and star-formation efficiency based on our numeri- 
cal simulations combined with the model reionization histo- 
ries suggest that the observed transmission and luminosity 
function are compatible with reionization at 2 ~ 6. 



7.2 Comparison of the escape fraction in 
numerical and semi analytic models 

Our immerical simulations are only able to resolve galax- 
ies down to a halo mass of M = 2 x IO^A/q, while in 
our fiducial reionization model the bulk of ionizing pho- 
tons prior to reionization come from lower mass galaxies 
(M > 10* M0). We therefore also use our semi-analytic his- 
tory model to investigate the implications of excluding these 
lower mass galaxies, as well as galaxies with halo mass below 
10^" M© (see § 2)). We adopt values for the escape fraction 
of /esc = 0.01 and f*^^ = 0.02 respectively, which are con- 
sistent with our numerical results in the previous sections. 
The results of these models are plotted in Figure [§] (thin 
and dashed lines). These histories are consistent with reion- 
ization at z ~ 6, which directly constrains the the result 
that suggests that the observed transmission and luminosity 
function are compatible with reionization at z ~ 6. However, 
these models exceed the observed ionization rate at z < 6, 
owing to the rapid increase in the collapsed fraction. 

In addition to the properties described in Figure |9l 
reionization histories are also constrained by the optical 
depth to Thomson scattering of CMB photons. The fiducial 
model in Figure[9]has an optical depth of Tos = 0.063, while 
models with Mesc = W^Mq and Mosc = 10^'^Mq have t^s = 
0.045, Tos ~ 0.041. These values are smaller than the value 
determined froin the WMAP satellite, Tea = 0.084 ± 0.016 
(|Komatsu et al.ll2008l ). 



7.3 Reionization histories with Massive Stars 

Before concluding, we consider the implications of a pop- 
ulation of massive stars at high redshift. Models which 
include massive population-Ill stars partially reionize the 
IGM at high redshift an d result in larger value s of Tes - For 
example, in th e work of IChoudhurv fc Ferraral l|2005r ) and 
IWvithe fc CenI (|2007l ) , models in which the reionization his- 
tory shows an extended plateau of mostly ionized IGM in 
the range 6 < z < 10 are able to produce values of Tes within 
the preferred range plus completion of reionization at z ~ 6. 
This early reionization does not however make our predicted 
values of the escape fraction inconsistent with reionization 
at z ~ 6. 

The heavy dashed line in Figure|9]represents a model for 
the reionization history which includes population-Ill stars 
at high redshift. Here we assume that the value of f*^^ is not 
sensitive to redshift, but that 10 times the number of ionizing 
photons per baryon are emitted prior to a transition redshift 
Ztran. Note that this modelling is intended to be illustrative 
only and does not include more advanced t reatments of the 
transi ti on redshift, such as those studied bv lSchneider et al.l 
(l2006t ). IScannapieco et all l|2003l ) and Wyithe fc Cen (2007). 
The example shown illustrates the effect of a high redshift 
population of massive stars and corresponds to an assumed 
transition redshift of Ztran = 8. Finally, the addition of 
population-Ill stars increases the value of Tos. For the in- 
clusion of population-Ill stars in our fiducial model, where 
Ztran ~ 8, we fiud Tos ~ 0.095, which is consistent with the 
resuhs from the WMAP satellite (tcs = 0.084 ± 0.016). 
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8 CONCLUSIONS 

In this paper we have employed an A'^-body simulation 
to describe the density and peculiar velocity fields in a 
(65.6 Mpc/;i)-' volume of the IGM aX z = 5.5, 5.7, 6.0. The 
simulation also identified the locations and masses of virial- 
ized halos, with the minimum halo mass accurately resolved 
by the simulation being 2 x 10^ Mq. We associate the identi- 
fied halos with galaxies and assign to them a UV luminosity 
using models for the S FR and the SEP. By fitting to the ob- 
served UV galaxy LF (JBouwens et al.ll2007l ) we constrained 
the best fit free parameters to be /* =0.11 and ^dc =0.16, 
where /* is the fraction of baryons which participate in star 
formation and ioc is the duty-cycle of star-bursting galaxies. 
We then constructed an ionizing background which reflects 
the spatial distribution of UV luminous galaxies. Our model 
made the assumption of a universal MFP at each redshift. 

Given the density, velocity and ionizing radiation flelds 
we extracted an ensemble of Lya absorption spectra along 
3072 sight-lines. From this sample we computed the mean 
effective transmission as a function of the fraction of ionizing 
photons produced that escape their host galaxy (/esc). In our 
simulations /esc effectively adjusts the intensity of the ioniz- 
ing background. The value of fosc was determined so that the 
mean effective transmission predicted by the model matched 
the observed value. Assuming halos above the mass resolu- 
tion of our simulation to host galaxies, we find values for 
the escape fraction at 2; ~ 5.5 — 6 of /esc -^ 10 — 2 5% assum- 
ing v alues for the MFP which are suggested by (jFan et al.l 
120061 ). In models where we assume the formation of galaxies 
to be suppressed in halos with masses below 10^" M0 (pos- 
sibly due to radiation feedback following reionization) , we 
find the escape fraction to be in the range /esc ~ 20 — 40%. 
Our res ults for the escape fraction are consistent with those 
found in lBolton fc Haehnclt (^2007) who constrain the escape 
fraction to be /eac ~ 20 — 30%, using a similar approach to 
the one adopted in this paper. This is true both for scenar- 
ios in which the minimum halo mass is set by the simulation 
resolution (2 x 10^ Mq) and for the larger value of 10^° Mq. 

We also considered the scenario suggested by 
iGnedin et al.l (|2007l ) in which the escape fraction of ion- 
izing photons is negligible in halos with masses below ~ 
5 X IO^^Mq. In this case we find the escape fraction is re- 
quired to be /esc ~ 46 — 86% in order to reproduce the ob- 
served optical depth of the IGM. This is inc onsistent with 
theore tical modelling of the escape fraction bv lGnedin et al.l 
(|2007|) who suggest /esc ~ 1 — 3%. Using the numerical re- 
sults to calibrate an analytic relation between the escape 
fraction and minimum galaxy halo mass we extrapolate our 
results to a mass {M ~ 10* Mq ) corresponding to the hy- 
drogen cooling threshold. In this case we find /esc ~ 5— 10%, 
consistent with observed estimates at lower redshift. 

Since the SFR is estimated at 1350A, while /esc is es- 
timated at 912^4, the escape fraction is degenerate with the 
adopted galaxy SED. The value of the Lyman-break factor 
(LB), which is the intrinsic flux decrement measured across 
the Lyman break, is particularly sensitive to to the age of 
the stellar population. Assuming that the normalization of 
the adopted SED at wavelengths above the Lyman-limit, 
and the slope of the SED at wavelengths below the Lyman- 
limit are both correct, we can express this degeneracy with 
respect to LB as /esc = ./^c"^ {LB/3). In this expression 



/csc"^ is the escape fraction which we have constrained in 
this paper, and /esc the value obtained for an alternate LB. 

We use a semi-analytic model that describes the evo- 
lution of the ionized state of the post-reionization IGM to 
study whether the escape fraction implied by our numer- 
ical modelling is consistent with reionization a,t z > 6. 
We find that the measured ionization rates imply a his- 
tory that is consistent with reionization at z ~ 6 for a 
range of these models, including those in which the min- 
imum halo mass of galaxies contributing to the ionizing 
emissivity is M ~ 10^,10^ and 10^" solar masses respec- 
tively. However, we find that models with minimum masses 
above 10^ Mq over-estimate the ionization rate in the post- 
overlap IGM. These models imply that the Universe was 
reionized by low mass galaxies with an escape fraction of 
-5%. 

In summary, the transmission of Lya photons through 
the IGM at 2 ~ 5.5 — 6 implies an ionization rate that re- 
quires the escape fraction of ionizing photons to be in excess 
of 5%. Our results for the escape fraction and star-formation 
efficiency based on our numerical simulations combined with 
the model reionization histories suggest that the observed 
transmission and luminosity function are compatible with 
reionization at z ~ 6. 
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